Polyfem Examples¶
This is a jupyter notebook. The "real" notebook can be found here.
Some imports for plotting
import plotly.offline as plotly
import plotly.graph_objs as go
#Necessary for the notebook
plotly.init_notebook_mode(connected=True)
algebra
import numpy as np
and finallypolyfempy
import polyfempy as pf
Plotting utility¶
This code is not particularly interesting.
It converts a tet-mesh (4 indices) into a face-based tri mesh and uses plotly Mesh3d to plot it.
def plot(vertices, connectivity, function):
x = vertices[:,0]
y = vertices[:,1]
if vertices.shape[1] == 3:
z = vertices[:,2]
else:
z = np.zeros(x.shape, dtype=x.dtype)
if connectivity.shape[1] == 3:
f = connectivity
else:
#Convert a tet-mesh into face based triangles.
f = np.ndarray([len(connectivity)*4, 3], dtype=np.int64)
for i in range(len(connectivity)):
f[i*4+0] = np.array([connectivity[i][1], connectivity[i][0], connectivity[i][2]])
f[i*4+1] = np.array([connectivity[i][0], connectivity[i][1], connectivity[i][3]])
f[i*4+2] = np.array([connectivity[i][1], connectivity[i][2], connectivity[i][3]])
f[i*4+3] = np.array([connectivity[i][2], connectivity[i][0], connectivity[i][3]])
mesh = go.Mesh3d(x=x, y=y, z=z,
i=f[:,0], j=f[:,1], k=f[:,2],
intensity=function, flatshading=function is not None,
colorscale='Viridis',
contour=dict(show=True, color='#fff'),
hoverinfo="all")
layout = go.Layout(scene=go.layout.Scene(aspectmode='data'))
fig = go.Figure(data=[mesh], layout=layout)
plotly.iplot(fig)
Set the mesh path
mesh_path = "plane_hole.obj"
create settings
settings = pf.Settings()
pick linear $P_1$ elements. If the mesh would be a quad it would be $Q_1$
settings.discr_order = 1
normalize the mesh to be in [0,1]^2
settings.normalize_mesh = True
and choose Young's modulus and poisson ratio
settings.set_material_params("E", 210000)
settings.set_material_params("nu", 0.3)
We are use a linear material model
settings.tensor_formulation = pf.TensorFormulations.LinearElasticity
Next we setup the problem
problem = pf.GenericTensor()
sideset 1 has zero displacement in $x$
problem.add_dirichlet_value(1, [0, 0], [True, False])
sideset 4 has zero displacement in $y$
problem.add_dirichlet_value(4, [0, 0], [False, True])
sideset 3 has a force (Neumann) of [100, 0] applied
problem.add_neumann_value(3, [100, 0])
fianally set the problem
settings.set_problem(problem)
Solve!
solver = pf.Solver()
solver.settings(str(settings))
solver.load_mesh(mesh_path)
solver.solve()
Get the solution
[pts, tets, disp] = solver.get_sampled_solution()
diplace the mesh
vertices = pts + disp
and get the stresses
mises, _ = solver.get_sampled_mises_avg()
finally plot with the above code
plot(vertices, tets, mises)
Note that we used get_sampled_mises_avg to get the Von Mises stresses because they depend on the Jacobian of the displacement which, becase we use $P_1$ elements, is piece-wise constant. To avoid that effect in get_sampled_mises_avg the mises are averaged per vertex weighted by the area of the triangles. If you want the "real" mises just call
mises = solver.get_sampled_mises()
plot(vertices, tets, mises)
Torsion¶
Non-linear example. We want to torque a 3D bar around the $z$ direction.
The example is really similar as the one just above.
Sets the mesh and create the settings. As before
mesh_path = "square_beam_h.HYBRID"
settings = pf.Settings()
It is an hex-mesh so we are using $Q_1$
settings.discr_order = 1
In this case the mesh has already the correct size.
settings.normalize_mesh = False
Choose Young's modulus and poisson ratio, as before
settings.set_material_params("E", 200)
settings.set_material_params("nu", 0.35)
Differently from before we want a non-linear material model: NeoHookean
settings.tensor_formulation = pf.TensorFormulations.NeoHookean
and we want to do 5 steps of incremental loading to avoid ambiguities in the rotation
settings.nl_solver_rhs_steps = 5
Now we setup problem with fixed sideset, rotating an number of tours
problem = pf.Torsion()
sideset 5 is fixed
problem.fixed_boundary = 5
sideset 6 rotates
problem.turning_boundary = 6
around the $z$-axis (2)
problem.axis_coordiante = 2
by half a tour
problem.n_turns = 0.5
Now we choose a coarse visualization mesh
settings.vismesh_rel_area = 0.001
and set the problem and solve
settings.set_problem(problem)
#solving!
solver = pf.Solver()
solver.settings(str(settings))
solver.load_mesh(mesh_path)
solver.solve()
takes approx 1 min because it is a complicated non-linear problem!
Get solution and stesses as before
[pts, tets, disp] = solver.get_sampled_solution()
vertices = pts + disp
mises, _ = solver.get_sampled_mises_avg()
and plot the 3d result!
plot(vertices, tets, mises)
</div>